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REMARKS 

In the Final Office Action dated 29 December 2005, claims 1-16 were pending; claims 
1-11 and 13-16 were rejected; and claim 12 was objected to but considered to include 
patentable subject matter. 

By way of this Amendment, claim 16 is amended to better represent the present 
invention. No new matter was added by way of this amendment. 

In the office action claims 1, 3-11, 13-14, and 16 were rejected under 35 USC § 103 (a) 
as being unpatentable over Trampert in view of Zhu. Further, claims 1, 3-8, 10, 13-14, and 16 
were rejected under 35 USC § 103 (a) as being unpatentable over Trampert in view of 
Silawongsawat. 

The present claim 1 includes "estimating a P-SV propagator from said recorded seismic 
wavefield". It is submitted that this step goes beyond the mere knowledge of a matrix 
representation of the P-SV propagator. This step requires an estimation of the P-SV 
propagator from the recorded seismic wavefield. 

In the following paragraphs, it will be shown that neither Trampert nor Zhu nor 
Silawongsawat teach a way of estimating a P-SV propagator from a recorded seismic 
wavefield. 

Regarding first Trampert, Trampert makes an explicit statement that the method 
presented is only applicable to the SH propagator matrix. With regard to the method used by 
Trampert, it is stated that "the spectral ratio of the total surface record over borehole record 
gives a measurement of the propagator" (see Trampert, p. 291, 2 nd col, 2 nd paragraph). The 
approach of Trampert is only valid for SH waves as such wave types do not undergo 
conversion at the surface. 

Regarding secondly Zhu, it has to be noted that Zhu is not concerned with an estimate 
of the P-SV propagator. Zhu takes the Thompson-Haskell expression of the P-SV propagator 
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(equation 20). The Thompson-Haskell expression of the P-SV propagator is equivalent to the 
expression given by Aki and others. Zhu shows that the P-SV converges to a static solution 
if the frequencies go to zero. Zhu shows further that static deformation from earthquakes can 
be derived using the P-SV propagator. However, the important point is that Zhu does not 
show how to estimate the S-PV propagator from a recorded seismic wavefield. Instead, Zhu 
uses a predefined earth model referred to as "PREM" (see page 5, 2 nd col, bottom paragraph) 
to define the parameters of the P-SV propagator. 

When finally considering Silawongsawat, it is again clear that as far as the formulation 
of the P-SV propagator is concerned, Silawongsawat recites the known expressions based on 
Backus, Thompson-Haskell, Aki and others. For the present discussion, it should be noted 
that Silawongsawat, similar to Zhu, teaches methods for numerical simulation. In these 
methods, the P-SV propagator is not estimated from recorded seismic data, but derived solely 
from the parameters which define the earth model used in the simulation. 

In the rejection of claims 1, 3-1 1, 13-14, and 16 over Trampert in view of Zhu or in 
view of Silawongsawat, the Examiner neglected an important feature of independent claims 1 
and 16. More specifically, the Examiner appears to equate the prior knowledge of an 
mathematical expression for the P-SV propagator with the estimation of a P-SV propagator 
from recorded seismic data. It is, however, evident that knowledge of the P-SV propagator is 
not equal to a method of estimating it from real data. The Applicants have not denied that 
mathematical expressions of the P-SV operator are prior knowledge. In fact equations (9) - 
(14) of the specification represent such an expression. Those equations are taken from Aki et 
al. and serve as a starting point for the estimation of the P-SV operator from recorded data. 

In this context, Applicants would like to comment that being restricted to claim 12, 
which was deemed to include allowable subject matter, is too narrow. The equations (17)- 
(21) denote only an intermediate step in the estimation of the P-SV operator. Once the terms 
of the P-SV propagator is split in real an imaginary parts (as per equation (17) - (21)), it can 
be quickly realized that, for example, for near vertical incidents of the waves, the main 
components of the propagator can be approximately estimated from the recorded horizontal 
and vertical waves only. This variant of the present invention is explained in the last 
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paragraph of page 19 of the specification. The amended claim 16 broadens the scope of 
claim 12 to include all the methods described on page 19. 

Therefore, Applicants respectfully traverse all of the section 103 rejections in the Office 
Action because the cited references do not teach or suggest all of the claim limitations of 
independent claims 1 and 16. Specifically, the Examiner has not shown a reference that 
describes estimating a P-S V propagator from recorded seismic data. 

Moreover, Applicants respectfully submit that the Trampert, Zhu and Silawongsawat 
references do not contain any suggestion or motivation for combining the references. Whilst 
Trampert relates to identifying parameters relating to near-surface layers Zhu and 
Silawongsawat, both relate to modelling. As such, it would require impermissible hindsight 
to combine the references. 

Further, it appears worth noting that since the January 15, 2003 GB filing date of the 
present application, two papers by the co-inventors on the method of the present invention 
have been published in peer-reviewed journals. This underlines the fact that the scientific 
community has accepted the present invention as representing novel and original work. The 
publications are 

1. R. van Vossen, A. Curtis and J. Trampert 2005. Subsonic near-surface P velocity and 
low S velocity observations using propagator inversion. Geophysics, 70(4), pages R15-R23; 
and 

2. R. van Vossen, J. Trampert and A. Curtis 2004. Propagator and wave-equation 
inversion for near-receiver material properties. Geophys. J. Int., 157, 796-812. 

Courtesy copies of these documents are provided herewith. 



Hence, it is respectfully requested that the section 103 rejections of 
independent claims 1 and claim 16 be withdrawn. Additionally, it is respectfully requested 
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that the section 103 rejection of all claims depending from independent claims 1 and 16 also 
be withdrawn. 

In light of the above remarks, Applicants believe that the present application and claims 
1-16 as amended are in proper condition for allowance. Such allowance is earnestly 
requested. If the Examiner is contemplating any action other than allowance of all pending 
claims, the Examiner is urged to contact Applicant's representative, Jody Lynn DeStefanis, at 
(203) 431-5505. 

Applicants believe that no fee is due in connection with this response. However, in the 
event that a fee or refund is determined to be due, the Commissioner is hereby authorized to 
charge any underpayment or credit any overpayment to Deposit Account No 19-0615. 



Schlumberger Doll Research Center 
36 Old Quarry Road 
Ridgefield, CT 06877-4108 
Phone: (203)431-5505 
Fax: (203)431-5640 
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Subsonic near-surface P-velocity and low S-velocity 
observations using propagator inversion 
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ABSTRACT 

Detailed knowledge of near-surface P- and S-wave ve- 
locities is important for processing and interpreting multi- 
component land seismic data because (1) the entire wave- 
field passes through and is influenced by the near-surface 
soil conditions, (2) both source repeatability and receiver 
coupling also depend on these conditions, and (3) near- 
surface P- and S-wave velocities are required for wave- 
field decomposition and demultiple methods. However, 
it is often difficult to measure these velocities with con- 
ventional techniques because sensitivity to shallow-wave 
velocities is low and because of the presence of sharp 
velocity contrasts or gradients close to the earth's free 
surface. We demonstrate that these near-surface P- and 
S-wave velocities can be obtained using a propagator in- 
version. This approach requires data recorded by at least 
one multicomponent geophone at the surface and an ad- 
ditional multicomponent geophone at depth. The propa- 
gator between them then contains all information on the 



INTRODUCTION 

Strong near-surface velocity contrasts are often encoun- 
tered in land seismic surveys. Both P- and S-wave velocities 
may increase by nearly an order of magnitude at the interface 
defining the top of the bedrock, and P-velocities increase up to 
100% across the top depth of total water saturation (Stiimpel 
et al., 1984; Goforth and Hayward, 1992). 

Detailed knowledge of near-surface velocities is essential 
for engineering applications and groundwater and environ- 
mental projects (Ward, 1990). Furthermore, this knowledge is 
required to correctly process and interpret (multicomponent) 
land data. For instance, near-surface soil conditions have a 
significant influence on source wavelet and radiation patterns 
(Kahler and Meissner, 1983; Aritman, 2001). Also, wavefield 



medium parameters governing wave propagation between 
the geophones at the surface and at depth. Hence, invert- 
ing the propagator gives local estimates for these parame- 
ters. This technique has been applied to data acquired in 
Zeist, the Netherlands. The near-surface sediments at this 
site are unconsolidated sands with a thin vegetation soil 
on top, and the sediments considered are located above 
the groundwater table. A buried geophone was positioned 
1.05 m beneath receivers on the surface. Propagator in- 
version yielded low near-surface velocities, namely, 270 ± 
15 m/s for the compressional-wave velocity, which is well 
below the sound velocity in air, and 150 ± 9 m/s for 
the shear velocity. Existing methods designed for imag- 
ing deeper structures cannot resolve these shallow material 
properties. Furthermore, velocities usually increase rapidly 
with depth close to the earth's surface because of increas- 
ing confining pressure. We suspect that for this reason, sub- 
sonic near-surface P-wave velocities are not commonly ob- 
served. 



decomposition, which enables independent interpretation of 
up- and downgoing P- and S-waves, requires the free-surface 
reflectivity to be known accurately (e.g., Dankbaar, 1985; 
Wapenaar et al., 1990; Robertsson and Curtis, 2002). Wave- 
field decomposition is a prerequisite for demultiple methods 
(Verschuur et al., 1992), which are especially important in me- 
dia with a near-surface low-velocity layer that may act as a 
wave guide in which energy may propagate over long distances 
with little loss from geometric spreading. This could mask re- 
flections from a deeper target (Hunter et al., 1984; Roberts- 
son et al., 1996). Demultiple methods remove these guided 
waves. 

While shallow material properties are especially important 
for processing and interpretation of multicomponent seismic 
data, near-surface wave velocities usually cannot be resolved 
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with an acquisition geometry designed for imaging deeper 
structure. Detailed information can, however, be obtained 
with shallow, high-resolution reflection and refraction exper- 
iments (Doornenbal and Helbig, 1983; Hunter et al., 1984; 
Steeples and Miller, 1990). These techniques use arrays of 
closely spaced geophones and high frequencies to obtain de- 
tailed images of the shallow subsurface. 

Subsonic P-wave velocities have been observed with these 
shallow, high-resolution seismic experiments by analyzing 
moveout velocities close to the source (Birkelo et al., 1987; 
Bachrach and Nur, 1998; Bachrach et al., 1998; Baker et al., 
1999), whereas they are not commonly observed using con- 
ventional seismic techniques. This is a consequence of the dif- 
ferent depth sensitivities of these methods, combined with a 
near-surface velocity gradient caused by increasing confining 
pressure. A drawback of estimating near-source moveout ve- 
locities is that the complexity of the near-source field requires 
careful processing and interpretation of these types of data to 
avoid misinterpretation of recorded events (Michaels, 2002). 

Recently, Curtis and Robertsson (2002) introduced a tech- 
nique for estimating local near-surface velocities using a 3D 
geophone configuration. Geophones are not only deployed at 
the surface but also at shallow depths to enhance imaging of 
the near-surface without having to perform an additional high- 
resolution experiment. With the proposed 3D geophone con- 
figuration, spatial wavefield derivatives can be approximated, 
allowing inversion of the wave equation for near-surface 
P- and S-wave velocities (Curtis and Robertsson, 2002). An 
advantage of this method is its applicability to the complete 
wavefield; a drawback is its sensitivity to deployment-related 
errors (Muijs et al., 2002). 

We present results from a technique referred to as propaga- 
tor inversion (PI) (Trampert et al., 1993; Van Vossen et al., 
2004). This technique also uses a 3D geophone configura- 
tion to determine near-surface P- and S-wave velocities, but 
it avoids explicit computation of spatial wavefield derivatives 
and is therefore less sensitive to deployment-related errors. 
Moreover, it does not require measurement and interpreta- 
tion of moveout velocities in the near-source region, and it can 
be incorporated in a seismic survey for imaging deeper struc- 
ture without having to perform an additional high-resolution 
experiment. 

PROPAGATOR ESTIMATION FROM DATA 

Propagator matrices were introduced in seismology by 
Thomson (1950) and Haskell (1953) and generalized by 
Gilbert and Backus (1965). These matrices describe the propa- 
gation of plane waves through a horizontally layered medium. 
Throughout this paper, the free surface is used as a ref- 
erence level. The propagator can then be interpreted as a 
wavefield-extrapolation filter. Application of the propagator 
to the recorded wavefield at the free surface gives the wave- 
field at depth Az. 

Trampert et al. (1993) introduced SH PI to obtain the SH- 
wave velocity structure and the quality factor in a borehole. 
This propagator can be obtained from the recorded data by 
taking the spectral ratio of a downhole record over a surface 
record. It is completely determined by the medium parameters 
governing wavefield propagation between these two records. 
Recently, Van Vossen et al. (2004) formulated propagator es- 
timation for the elastic P-SV case. We briefly review this con- 



cept before we discuss the inversion scheme for near-surface 
material parameters. 

In an isotropic medium, the propagator naturally decom- 
poses into SH and coupled P-SV waves. The anelastic SH 
case is fully treated by Trampert et al. (1993) We only review 
the elastic P-SV case here. Denote the inline particle-velocity 
component by v\ and the vertical component by u 3 . The full 
propagator is a 4 x 4 matrix, and the boundary conditions 
state that the free surface is stress free, so that the wavefield 
at depth Az is related to the wavefield recorded at the free 
surface (z = 0) by: 



v\(co, x y Az) 




V\((O i x, 0) 
V 3 ((D,X,0) 



(1) 



For an elastic, homogeneous medium, with P velocity a and S 
velocity fi, the propagator coefficients read in the time domain 
(Aki and Richards, 2002; Van Vossen et al., 2004) as 

P n = P 2 P 2 G? + [(1 - 2/*V)/2]G? , (2) 

P 33 = [(1 - 20 V)/2]Gf + p 2 p 2 Gl (3) 

P 13 = [p(l - 2fi 2 p 2 )/(2q P )} G$ - p 2 pq s G s 2 , (4) 

P 31 = fpqrG* - [ P {\ - 2(i 2 p 2 )/(2q s )] Gf , (5) 

where 

Gf = 8(t + q P Az) + 8{t - q P Az), (6) 

G 2 P = 8(t + q P Az) - 8{t - qpAz), (7) 

Gf = 8(t + q s Az) + 8(t - q s Az), (8) 

G s 2 = 8(t + q s Az) - 8(t - q s Az). (9) 

The horizontal slowness is denoted by p, and the vertical slow- 
nesses q P and q s are 



qp = {a-*-p 2 y\ 

qs = {p- 2 -p 2 ) 112 . 



(10) 
(11) 



These theoretical expressions show that P n and P33 are sym- 
metric around t = 0, whereas P13 and P31 are antisymmetric 
around t = 0. Thus, in the frequency domain, P\\ and P33 are 
entirely real, and Pn and P 3J are purely imaginary. As a result, 
we can directly estimate the components of Pij(co) by equating 
real and imaginary parts in equation 1. In the following, we de- 
note the propagator coefficients estimated from the data with 
P(«t>). The explicit expressions for estimating P(o;) are 

P n =i®[v 3 ((o,0)p[v l (co,Az)] 

+ S[v 3 (a>, 0)p[ Vl (co, Az)]}/D(a>), (12) 
P33 = {9iMw,0)]!R[v3(w, Az)] 

+ SMa>, 0)]3[t;3(«, Az)]}/D(co), (13) 
A3 = iWlvxico, (^[ujK Az)] 

-SM^O)]^,^, Az)]}/D(w), (14) 
h\ =/{9f[u3(w,0)]3[t> 3 (w, Az)] 

- S[v 3 (o>, 0)]m[v 3 (a>, Az)]}/D(a>), (15) 
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with the denominator D(co) given by 

D(a>) = m[v 3 ((o, 0)]9t[ui(a>, 0)] + %[v 3 (<o, 0)p[vi(w, 0)]. 

(16) 

In these equations, 9t[ui(<a, Az)] is the real part of v\(co, Az), 
3[ui(w, Az)] is the imaginary part, and * denotes the imaginary 
unit Note that the symmetry properties used to obtain 
explicit expressions for the propagator filters break down in 
the viscoelastic case. Then, only the SH propagator can be di- 
rectly obtained from the data (Trampert et ai., 1993). 

Van Vossen et al. (2004) computed P with a stabilized spec- 
tral division using the so-called water-level method (Helm- 
berger and Wiggins, 1971). Stabilization is required because of 
a limited bandwidth of D(co) f and interfering waves may intro- 
duce internal notches in D(a>). However, a problem associated 
with this method is that the amount of stabilization can influ- 
ence the estimates for the propagator (Ammon, 1991), which 
may affect the inversion results as well (Van Vossen et al., 
2004). To avoid these problems, we decided instead to imple- 
ment the spectral divisions in the time domain using a Wiener 
deconvolution scheme. Either symmetric (P u and P33) or an- 
tisymmetric filters (P13 and P 3i ) around t = 0 are constructed 
with N independent coefficients. Details on the implemen- 
tation of the acausal Wiener deconvolution can be found in 
Appendix A. 

PROPAGATOR INVERSION (PI) 

In the previous section, we showed that P can be obtained 
from data recorded by one surface geophone and one at depth. 
In this section, we outline the inverse procedure for estimating 
the near-surface velocities a and f$ from P. 

A flow diagram for the inverse problem is shown in Fig- 
ure 1. The 3D receiver configuration used for the Zeist field 
experiment in the Netherlands, conducted to test propagator 
estimation and inversion, is shown in Figure 2. We discuss this 
experiment in detail in the next section. The configuration has 
multicomponent geophones positioned in a cross shape at the 
surface. Two geophones are buried at the center of the receiv- 
ing group. Strictly speaking, propagator estimation requires 
only one multicomponent geophone positioned at the surface 
and a second geophone at depth. Then, the horizontal slow- 
ness also needs to be constrained by the inverse procedure. 



Data selection 



Data-esli mated 
propagator 



P a>P 



Theoretical 
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Frequency filtering 



I 



Waveform fitting 



Grid search 



However, we do not have to incorporate horizontal slowness 
in the inverse procedure when a short array of geophones is 
deployed at the surface, since it can be measured directly. 

PI consists of the following steps. First, a data window is se- 
lected to isolate an arrival. The window is tapered at its edges 
by a cosine taper. Second, the selected data are used as input 
for propagator estimation. The horizontal slowness p of the 
dominant arrival in this time window can be determined using 
the array of geophones in the inline direction x. This is accom- 
plished by estimating the time shifts for which the stack power 
is optimized. The remaining unknown parameters in the theo- 
retical propagator for a homogeneous isotropic medium are 
a and ft (equations 2-11). Values for a and are selected 
using a grid-search technique. A physical bound is imposed 
such that p< a/\/2, i.e., the Poisson's ratio has to be positive. 
Given values for a, /?, and /?, the theoretical propagator can be 
evaluated. Before comparing the waveforms of the theoretical 
propagator P to the data-estimated propagator P, frequency 
filtering is necessary since P is band-limited, whereas P has an 
infinite bandwidth. After bandwidth equalization, the propa- 
gator waveforms can be compared to each other. The L2 norm 
was used as the objective function: 



( NAt \ 
SiJ = \ E [Pij(t)-Pij(t,<*,frp)] 2 \ 



1/2 



(17) 



with i,j =1,3. The objective function for the joint inverse of 
all propagator coefficients is given by the sum of all individual 
misfit functions, 

E tot = E n + E n + £31 + £33. (18) 
Estimates for a and f$ can be obtained by minimizing E tot . 

APPLICATION ON ZEIST FIELD DATA 

We illustrate PI on a field data set acquired in Zeist, the 
Netherlands. On this site, the near-surface material mainly 
consists of unconsolidated sands, with a thin layer of vegetal 
soil on top. 
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Figure 1. Propagator inversion scheme for a and p. 



Figure 2. 3D receiver configuration for the Zeist field exper- 
iment, (a) top view and (b) front view. All geophones were 
multicomponent geophones, and the source positions were lo- 
cated on the x-axis. 
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Data acquisition 

A walkaway noise test was performed with 3C 4.5-Hz geo- 
phones at offsets between 0.75 and 84 m, with 0.75-m geo- 
phone spacing. These data (Figure 3) show that ground roll 
and guided waves are dominant in the recordings. In addi- 
tion, measurements were made with a dense 3D 3C receiv- 
ing configuration that will be used for PI (see Figure 2). Geo- 
phones were positioned in a cross shape at the surface. In the 
jc-direction, receivers were located at 0.50, 1.00, and 2.00 m 
distance to the center of the configuration; in the ^-direction, 
the distances were 0.15, 0.50, 1.00, and 2.00 m. The buried 
geophones were positioned at 0.45 and 1.05 m depth, respec- 
tively. Geophones can be buried efficiently in unconsolidated 
sediments using a hand ground drill. This approach minimizes 
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Figure 3. Walkaway noise survey for (a) v x and (b) v z . The panels are displayed 
with trace normalization. 



Table 1. Horizontal slowness estimates for each shot 
position. Data are selected in time windows with t between t\ 
and t 2 . 



Shot 


Offset 


h 


h 


P 




number 


(m) 


(s) 


W 


(ms/m) 


(ms/m) 


1 


85 


0.19 


0.26 


2.10 


0.11 


2 


80 


0.18 


0.25 


2.12 


0.07 


3 


75 


0.17 


0.24 


2.23 


0.03 


4 


70 


0.16 


0.23 


2.18 


0.06 


5 


65 


0.15 


0.22 


2.20 


0.06 


6 


60 


0.14 


0.21 


2.15 


0.06 


7 


55 


0.13 


0.20 


2.12 


0.04 


8 


50 


0.11 


0.18 


2.11 


0.03 


9 


45 


0.09 


0.17 


2.23 


0.05 


10 


40 


0.08 


0.16 


2.19 


0.05 


11 


35 


0.07 


0.14 


2.21 


0.03 


Mean 








2.17 


0.05 




0.1 0.12 
Time (s) 



0.14 



Figure 4. (a) Inline- and (b) vertical-component recordings 
with the source located at 35 m distance to the center of the 
receiver group. 



the medium perturbations caused by burial of geophones. For 
the dense 3D 3C group, data were acquired using 11 differ- 
ent source positions located between 35- and 85-m offset. The 
shot spacing was 5.0 m, and all shot points were located on 
the *-axis. The experiment was repeated four times for every 
shot position. During the whole experiment, a weight drop was 
used; a steel ball of about 37 kg was dropped from approx- 
imately 3.5 m height on a steel plate resting on the ground. 
The recording instrument was a Bison Spectra with 48 chan- 
nels, and the time-sampling interval was 0.25 ms. 

Data selection and estimation of horizontal slowness 

In theory, the propagator method is valid for a single slow- 
ness, i.e., isolating an arrival as uncontaminated as possible by 
other arrivals. On the other hand, consid- 
ering only very short time windows results 
in a poor signal-to-noise ratio. Synthetic ex- 
periments demonstrate that PI is insensi- 
tive to small slowness variations in the se- 
lected data (Van Vossen et al., 2004). There- 
fore, we selected data windows including all 
events arriving before the ground roll. For 
each shot position, time windows that con- 
tain the selected data are listed in Table 1. 
A cosine taper with 0.01 s length was ap- 
plied to both edges of the window. An ex- 
ample of data selected for slowness estima- 
tion is shown in Figure 4. A bandpass fil- 
ter with cut-off frequencies between 40 and 
140 Hz was applied to these recordings. The 
events shown may be interpreted as trapped waves above the 
groundwater table. Because there are significant differences 
in the recorded amplitudes on the inline component, we de- 
cided to estimate p using only the vertical component of the 
recorded particle velocity. These slowness estimates are given 
for each shot position in Table 1. The differences between the 
estimated horizontal slownesses for the different source posi- 
tions are small. 

Propagator estimation and inversion 

Contrary to horizontal-slowness estimation, no frequency 
filtering was applied to the selected data prior to propagator 
estimation. An example of the data used for propagator esti- 
mation is shown in Figure 5. It shows that recordings rapidly 
change with depth, especially on the vertical component. At 
0.45 m depth, high frequencies are strongly attenuated com- 
pared to the recordings obtained at the free surface and at 
1.05 m depth. On the other hand, the low-frequency content 
of the signal decays with depth. This characteristic behavior is 
caused by interference between the free-surface incident wave 
and its reflected and converted waves. 

Interpretation of the data recorded on the horizontal com- 
ponent is difficult. The two surface geophones show significant 
amplitude differences in the high-frequency part of the spec- 
trum. For lower frequencies, on the other hand, there is ex- 
cellent agreement between these recordings. This observation 
could indicate coupling differences between these two surface 
geophones (Krohn, 1984). Another interesting observation is 
that the amount of energy recorded in the 50- to 60-Hz fre- 
quency band is small, although we did not apply a notch filter 
to the data. 
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We demonstrate PI first using surface geophone 1 and the 
buried geophone at 1.05 m depth. Then, we consider the data 
recorded by surface geophone 2. Given the frequency content 
of the signal and the time-sampling interval of 0.25 ms, accu- 
rate velocity estimation with the geophone buried at 0.45 m 
depth is, in our opinion, not feasible. 

Figure 6 shows the data-estimated propagator for each shot 
position. The frequency passband is between 40 and 140 Hz. 
The theoretical propagator P is shown with a and p for which 
E tot is minimized. There is good agreement between P and P 
for most individual shots. Since the changes in the estimated 
horizontal slowness are small (Table 1), stacking of the data- 
estimated propagator components over all source positions is 
not in conflict with the single slowness assumption. This pro- 
cess enhances the signal-to-noise ratio of the data-estimated 
propagator. Figure 6 shows that an excellent fit is obtained 
between the averaged propagator components and the best- 
fitting theoretical propagator. 

The constraints offered by each individual propagator com- 
ponent are shown in Figure 7a-d. Misfit functions are shown 
for the stacked propagator, and the minima of the prestack 
propagators illustrate the uncertainty. No computations are 
performed with combinations for a and p for which the 
Poisson's ratio becomes negative. The misfit functions show 
that P n and P 3 i are dominantly sensitive to variations in 
whereas is more sensitive to variations in a. Pn contains 
information on both a and p. Because P33 is dominantly sen- 
sitive to variations in a, (ftp) 2 1 (equation 3). This is con- 
firmed by Table 1. Thus, for near- vertical incident waves, the 
PI is dominantly sensitive to phase differences rather than am- 
plitude effects as a result of interaction of the incident wave- 
field with the free surface. However, close to the critical angle 
for incident S-waves, the amplitude coefficient of Pn changes 
rapidly, which results in sensitivity for both ot and 
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Figure 5. Traces and amplitude spectra for v x (a, b) and for v z 
(c, d) for different depths, recorded at 35-m offset. The black 
and blue traces are the recordings acquired at the free surface 
by the geophones labeled 1 and 2 (Figure 2); the red traces 
were acquired at 0.45 m depth and the green traces at 1 .05 m 
depth. 



Figure 7e illustrates the joint inversion for all propagator 
coefficients. Both a and ft are well constrained. The minimum 
of the joint inversion is equal to the average of all minima of 
the misfit functions for each individual shot. We obtain the 
following estimates: a = 270 ± 15 m/s and p = 150 ± 9 m/s. The 
uncertainties given are the standard deviations of the variation 
of best estimates for the individual shots. 

So far we have only discussed the data-estimated prop- 
agator obtained using surface geophone 1. Because there 
are significant differences between the horizontal component 
recordings 1 and 2 (Figure 5), it is important to assess the con- 
sequences of these data differences on propagator estimation 
and inversion. This allows us to determine whether PI is robust 
in the presence of realistic data errors. Figure 8 shows that the 
match between P and P is not as good compared to the re- 
sults for surface geophone 1, although the fit is good for the 
first three shot points and for F 33 . This propagator coefficient 
is not significantly affected by data variations on the horizon- 
tal component: it relates the vertical component acquired at 
the free surface to the same component at depth (equation 1). 
The high-cut frequency was lowered to 100 Hz to reduce the 
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Figure 6. Individual and averaged (stacked) data-estimated 
propagators (solid) compared to the best fitting theoretical 
propagator (dashed). Shown are (a) P n , (b) P U) (c) Pu, and 
(d) P33. The data-estimated propagators are computed using 
surface geophone 1 and the buried geophone at 1.05 m depth. 
Frequency filters are applied with a passband between 40 and 
140 Hz. 
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effects of coupling errors to the velocity estimation. The mis- 
fit function for the joint inverse of all stacked propagators is 
shown in Figure 7f. The velocities corresponding to the mini- 
mum of E tot are a = 230 m/s and = 155 m/s. Averaging all 
individual shots gives a = 244 ± 20 m/s and ft = 152 ± 14 m/s. 
The estimates for agree well with the previously obtained 
velocity estimates, whereas a is less consistent. For shotpoints 
1, 2, and 3, which show a good match between P and P, we find 
that a > 260 m/s. 

Although geophones 1 and 2 were close together (0.30 m 
between them), the horizontal-component data recorded 
by these two geophones significantly differ for frequencies 
above 70 Hz. These amplitude differences are attributed to 
geophone-ground coupling. This refers to the accuracy with 
which a geophone measures the actual ground motion. It is 
especially relevant for horizontal-component recordings. A 
well-coupled horizontal geophone has a coupling-resonance 
frequency of 130 Hz, whereas poorly coupled horizontal geo- 
phones could have significantly lower (down to 30 Hz) reso- 
nance frequencies (Krohn 1984). For frequencies much lower 
than the coupling-resonance frequency, the geophone accu- 
rately follows the ground motion. 

The results of surface geophone 1 are not sensitive to 
changes in the high-cut frequency from 100 up to 140 Hz. In- 




Figure 7. Waveform misfit functions (a) £ u , (b) £13, (c) £31, 
(d) £33 for poststack P u , Pu, ^31, and P33 computed usingsur- 
face geophone 1. The combined constraints are shown in (e). 
The joint inverse results using surface geophone 2 are shown 
in (f). Contours are drawn for £;; = min(E t j)+c, with c = 0.02, 
0.05, and 0.10, and the indices 1 and j take the values 1 or 3. 
The open circle denotes the minimum of £,;, and the triangles 
indicate the positions of the minima of Eq for each individ- 
ual shot. No computations are performed in the nonphysical 
region p > a/y/2. 



creasing this frequency reduces the uncertainty in the velocity 
estimates. For surface geophone 2, on the other hand, increas- 
ing this frequency resulted in poor fits of the propagator co- 
efficients and large uncertainties in the estimated velocities. 
Therefore, we decided to use different frequency bands for 
the application of PI to data acquired by surface geophones 1 
and 2. The obtained results are barely influenced by the choice 
of the low-cut frequency. Lowering this frequency to 20 Hz 
yielded similar velocity estimates, although the data misfit be- 
tween P and P increased somewhat, resulting in larger un- 
certainties attached to these estimates. Therefore, we selected 
40 Hz as the low-cut frequency. 

Thus, the analysis indicates that best results are obtained 
with the data recorded by surface geophone 1: inversion re- 
sults are stable up to 140 Hz, a better data fit is obtained, and 
the data have a better resolving power for near-surface P- and 
S-wave velocities. Although the data quality of the recordings 
of surface geophone 2 is poorer, the estimate for the S-wave 
velocity is in agreement with the results obtained with surface 
geophone 1, and the difference between the obtained P-wave 
velocities is approximately 10% of the estimated value. This 
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Figure 8. Individual and averaged data-estimated propaga- 
tors (solid) compared to the best-fitting theoretical propagator 
(dashed). Shown are (a) Pu, (b) /> 13> (c) P 31 , and (d) P33. The 
data-estimated propagators are computed using surface geo- 
phone 2 and the buried geophone at 1.05 m depth. Frequency 
filters are applied with a passband between 40 and 100 Hz. 
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Figure 9. Near-source traces of walkaway noise survey for (a) 
V*, (b) v x with 0.005 s AGC window, (c) v z , and (d) v z with 
0.005 s AGC window. The airwave is indicated by the solid 
line, the groundwater refraction by the dashed line, and the 
arrows indicate the reflected P-wave at the groundwater table. 



indicates that the propagator inversion is robust in the pres- 
ence of measurement errors. 

DISCUSSION 

Low near-surface velocities are obtained with PI, namely, 
a = 270 ± 15 m/s and p = 150 ± 9 m/s using geophone 1. The 
Poisson's ratio a that corresponds to these velocities is 0.28 
with an uncertainty range between 0.18 and 0.34. Because the 
Poisson's ratio is sensitive to perturbations in the estimated 
velocities, it is difficult to make a sensible lithological inter- 
pretation. Despite this uncertainty, we may argue that the ob- 
served Poisson's ratio lies between the end-member models 
for dry, gas-saturated sands with a e [0.0 0.22] and water- 
saturated sands with a e [0.38 0.50] (Bourbie et al., 1987), 
which qualitatively makes sense because the considered sedi- 
ment was partially water-saturated. 

To test the results of PI, we analyzed the Zeist walka- 
way noise spread with recordings between 0- and 84-m offset. 
Near-offset sections of these data are shown in Figure 9. The 
receiver spacing is 0.75 m. No frequency filtering was applied 
to these multicomponent recordings. The airwave is clearly 
visible on both the inline and vertical components in the auto- 
matic gain control (AGC) plots. In the offset range between 0 
and 10 m, we do not observe coherent energy arriving before 
the airwave. Events with a higher moveout velocity arrive just 
after the airwave. This indicates that the near-surface veloc- 
ity is low, and that velocities increase with depth close to the 
free surface. Between 15- and 30-m offset, the refracted wave 
from the water table can be observed, and also the reflected 



wave from this interface can easily be identified on the verti- 
cal component. Thus, the near-offset section of the walkaway 
noise spread is qualitatively in agreement with the P-wave ve- 
locity obtained with PI for the very shallow near-surface. 

In comparision to the results obtained with PI, we found 
a higher P-wave velocity with a dispersion analysis of guided 
waves. This difference can be attributed to the different depth 
sensitivity of dispersion analysis. Propagator inversion is only 
sensitive to wave velocities between the free surface and the 
buried geophone, whereas dispersion analysis is sensitive to 
velocities in the entire layer above the water table. For this 
reason, these low P- and S-wave velocities are not commonly 
observed in seismic surveys with a deeper target. 

Although shallow material properties may influence the ac- 
tual wave propagation, we believe that these have most im- 
pact on the measured wavefield. Both energy transmitted into 
the ground and recordings of the wavefield by geophones de- 
pend on the near-surface soil conditions. The repeatability of 
the source mostly depends on the (an)elastic properties of the 
soil (Aritman, 2001). Both the amount of energy radiated into 
the ground and the radiation pattern are influenced by the 
near-source material properties. Lateral variations in near- 
surface material properties could lead to poor repeatability of 
the source, degrading the quality of the seismic section. In ad- 
dition, the measurements of the recording instruments differ 
from the actual ground motion. The so-called receiver cou- 
pling is influenced by the stiffness of the soil (Krohn, 1984). 
As a consequence, differences between recordings of adjacent 
receivers can exist as a result of coupling differences. Coupling 
errors especially affect the quality of converted-wave data. 
Thus, correct processing and interpretation of seismic data ac- 
quired in a land seismic survey require an understanding of 
both source and receiver coupling effects, and such an under- 
standing may be facilitated by the very near-surface velocity 
estimates obtained using PI. 

Another issue is the interaction of the wavefield with the 
free surface. Usually, we are only interested in the free-surface 
incident P- and S-waves, while receivers placed on land mea- 
sure the interaction of these incident wavefields with the free 
surface. In principle, we can obtain the free-surface incident 
P- and S-waves using wavefield decomposition. Wavefield 
decomposition requires as input the free-surface reflectivity, 
which depends on the P- and S-wave velocities just below the 
free surface. Because we measure the wavefield exactly at the 
free surface, we should not use effective or averaged medium 
parameters for wavefield decomposition. Thus, the very shal- 
low material properties obtained with PI may improve wave- 
field decomposition and demultiple methods. 

Although application of the propagator method would re- 
quire additional effort in acquiring data, we demonstrate that 
the technique provides additional information relevant for a 
land seismic survey. The propagator method can be incorpo- 
rated into a seismic survey without having to perform an addi- 
tional high-resolution experiment. 

Finally, it should be mentioned that propagator estimation 
assumes that the medium is elastic. Although attenuation can 
be significant in the weathered layer, it does not affect the ob- 
tained results, since the dominant wavelength of the analyzed 
signal is not much smaller than the distance between the sur- 
face and the buried geophone. 
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CONCLUSIONS 

The data-estimated propagator contains all information on 
the material parameters governing wave propagation between 
the free surface and the depth of a buried geophone. We 
applied PI on Zeist field data to determine the local near- 
surface velocities. This inversion yielded subsonic compres- 
sional wave velocity, a = 270 ± 15 m/s, and a low shear ve- 
locity, = 150 ± 9 m/s for the top meter. 

Although very shallow anomalies are considered to have 
a small impact on the wavefield propagation, these may sig- 
nificantly influence the wavefield recordings. Both the energy 
transmitted into the subsurface by the source and receiver 
coupling depend on very shallow material properties. Hence, 
lateral changes in material properties could lead to poor re- 
peatability of the source and receiver coupling differences. 
Also, corrections for the interaction of the wavefield with 
the free surface require these shallow wave velocities to be 
known. 
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APPENDIX A 

PROPAGATOR ESTIMATION 
USING WIENER DECONVOLUTION 

In this section, we demonstrate the procedure to obtain the 
propagator filters, which are either symmetric or antisymmet- 
ric about t = 0, using Wiener deconvolution. We closely follow 
Yilmaz (2001) in our derivation, with the exception that we es- 
timate acausal filters with symmetry conditions around t = 0. 

Suppose that f(t) and g(t) are continuous signals, and that 
h(t) is given by the deconvolution of f(t) by g(t): 

h(t) = f(t)*g(tr\ (A-l) 



or, equivalently, 



s(0 **(') = /('). 



(A-2) 



(A-4) 



where * is the convolution operator. The function h(t) rep- 
resents here the unknown propagator. Crosscorrelating equa- 
tion A-2 with g(t) gives 

R(t)*h(t) = q(t), (A-3) 

where /?(/) denotes the autocorrelation of g(t) and 

/oo 
/(/ + r)«(r)dr. 
OO 

Assume that we may approximate the propagator h(t) by a fil- 
ter with 2M + 1 independent coefficients. This reads, denoting 
the time series with a vector, 

h = [h-M /i-M-H . . . A-i Ao *i • • - h M -\ h M ] T , (A-5) 

where T is the transpose operator. For a correlation with a 
maximum correlation length of /V + 1, q has 2N + 1 coeffi- 



cients, and reads 

q = [q. N q-N+\ ■ • • <?-l <?0 q\ • . - QnY - (A-6) 

Then, we may recast equation A-3 in a discrete form: 

Rh = q. (A-7) 

The coefficients of the (2N + 1) x (2M + 1) autocorrelation 
matrix R are given by 



(A-8) 



where r k denotes the fcth lag of the autocorrelation of g, and 
N + 1 is the maximum correlation length. 

In order to take the symmetry conditions into account, we 
partition q, h, and R; 



-(::)■ 



q = 



and 



R = 



R" " R" + 



The partitioned vectors are given by: 

h- = [Ao/2 A-i A-2 ••■ A.a/F, 
[A 0 /2 h x h 2 A a/] 7 , 



q" = ko q~\ q-2 
q + = too q\ <?2 

and the submatrices of R read as 
/ r 0 n r 2 

n r 0 n 

R + ' + = R~ - — r 2 r l r 0 
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Note that an additional row and column are added in the sys- 
tem of equations because the coefficients qo and /i 0 /2 appear 
both in the positive and negative parts of the partitioned vec- 
tors. It can be verified that equation A-7 can be rewritten in 
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partitioned form as 




(A-17) 



For a symmetric filter, h" = h + , hence equation A-17 reduces 
to 

(R + ' + + R + )h+ = q ^ q t (A-18) 

whereas for an antisymmetric filter, h" = — h + . This gives 

(R + ' + - R + '") h+ = q+ ~ q . (A-19) 

These systems of equations can be solved for the indepen- 
dent filter coefficients using a damped least-squares solution. 
Prewhitening of the data is essential to avoid artifacts resulting 
from limited bandwidth. 



REFERENCES 

Aki, K., and P. G. Richards, 2002, Quantitative seismology, 2nd ed.: 

University Science Books. 
Ammon, C. J., 1991, The isolation of receiver effects from teleseismic 

P waveforms: Bulletin of the Seismological Society of America, 81, 

2504-2510. 

Aritman, B. C, 2001, Repeatability study of seismic source signatures: 

Geophysics, 66, 1811-1817. 
Bachrach, R., and A. M. Nur, 1998, High-resolution shallow-seismic 

experiments in sand, Part I: Water table, fluid flow, and saturation: 

Geophysics, 63, 1225-1233. 
Bachrach, R., J. Dvorkin, and A. M. Nur, 1998, High-resolution 

shallow-seismic experiments in sand, Part II: Velocities in shallow 

unconsolidated sand: Geophysics, 63, 1234-1240. 
Baker, G. S., D. W. Steeples, and C. Schmeissner, 1999, In-situ, high- 
resolution P-wave velocity measurements within 1 m of the earth's 

surface: Geophysics, 64, 323-325. 
Birkelo, B. A., D. W. Steeples, R. D. Miller, and M. Sophocleous, 

1987, Seismic reflection study of a shallow aquifer during a pumping 

test: Ground Water, 25, 703-709. 
Bourbie, T., O. Coussy, and B. Zinszner, 1987, Acoustics of porous 

media: Technip. 

Curtis, A., and J. O. A. Robertsson, 2002, Volumetric wavefield 
recording and wave equation inversion for near-surface material 
properties: Geophysics, 67, 1602-1611. 

Dankbaar, J. W. M., 1985, Separation of P- and S-waves: Geophysical 
Prospecting, 33, 970-986. 

Doornenbal, J. C., and K. Helbig, 1983, High-resolution seismics on 
a tidal flat in the Dutch Delta — Acquisition, processing and inter- 
pretation: First Break, 1, 9-20. 



Gilbert, F., and G. E. Backus, 1965, Propagator matrices in elastic 
wave and vibration problems: Geophysics, 31, 326-332. 

Goforth, T., and C. Hayward, 1992, Seismic reflection investigations 
of a bedrock surface buried under alluvium: Geophysics, 57, 1217- 
1227. 

Haskell, N., 1953, The dispersion of surface waves in multilayered me- 
dia: Bulletin of the Seismological Society of America, 43, 17-34. 

Helmberger, D., and R. A. Wiggins, 1971, Upper mantle structure of 
Midwestern United States: Journal of Geophysical Research, 76, 
3229-3245. 

Hunter, J. A., S. E. Pullan, R. A. Burns, R. M. Gagne, and R. L. 
Good, 1984, Shallow seismic reflection mapping of the overburden 
bedrock interface with the engineering seismograph — Some sim- 
ple techniques: Geophysics, 49, 1381-1385. 

Kahler, S., and R. Meissner, 1983, Radiation and receiver pattern 
of shear and compressional waves as a function of Poisson's ratio: 
Geophysical Prospecting, 31, 421-435. 

Krohn, C. E., 1984, Geophone ground coupling: Geophysics, 49, 722- 
731. 

Michaels, P., 2002, Identification of subsonic P-waves: Geophysics, 67, 
909-920. 

Muijs, R., K. Holliger, and J. O. A. Robertsson, 2002, Perturbation 
analysis of an explicit wavefield separation scheme for P- and S- 
waves: Geophysics, 67, 1972-1982. 

Robertsson, J. O. A., and A. Curtis, 2002, Wavefield separation using 
densely deployed three-component single-sensor groups in land 
surface-seismic recordings: Geophysics, 67, 1624-1633. 

Robertsson, J. O. A., K. Holliger, A. G. Green, A. Pugin, and R. D. 
Iaco, 1996, Effects of near-surface wave guides on shallow high- 
resolution seismic refraction and reflection data: Geophysical Re- 
search Letters, 23, 495^198. 

Steeples, D. W., and R. D. Miller, 1990, Seismic reflection methods 
applied to engineering, environmental and groundwater problems, 
in S. H. Ward, ed., Geotechnical and environmental geophysics I: 
SEG, 1-30. 

Stiimpel, H., S. Kahler, R. Meissner, and B. Milkereit, 1984, The use of 
seismic shear waves and compressional waves for lithological prob- 
lems of shallow sediments: Geophysical Prospecting, 32, 662-675. 

Thomson, W. T., 1950, Transmission of elastic waves through a strati- 
fied solid medium: Journal of Applied Physics, 21, 89-93. 

Trampert, J., M. Cara, and M. Frogneux, 1993, SH propagator matrix 
ana Q s estimates from borehole- and surface-recorded earthquake 
data: Geophysical Journal International, 112, 290-299. 

Van Vossen, R., J. Trampert, and A. Curtis, 2004, Propagator and 
wave-equation inversion for near-receiver material properties: 
Geophysical Journal International, 157, 796-812. 

Verschuur, D. J., A. J. Berkhout, and C. P. A. Wapenaar, 1992, Adap- 
tive surface-related multiple elimination: Geophysics, 57, 1166- 
1177. 

Wapenaar, C. P. A., P. Herrmann, D. J. Verschuur, and A. J. 
Berkhout, 1990, Decomposition of multicomponent seismic data 
into primary P- and S-wave responses: Geophysical Prospecting, 38, 
633-661. 

Ward, S. H., 1990, Geotechnical and environmental geophysics, 
vol. I— III: SEG. 

Yilma2, O., 2001, Seismic data analysis I — Processing, inversion and 
interpretation of seismic data: SEG. 



2004 16:48 Geophysical Journal International 
. Geophys. J. Int. (2004) 157, 796-812 



gj 12249 



doi: 1 0. 1 1 1 1 /j. 1 365-246X.2004.02249.X 



Propagator and wave-equation inversion for near-receiver 
material properties 

R. van Vossen, 1 J. Trampert 1 and A. Curtis 2,3 

1 Faculty of Earth Sciences, Utrecht University, Utrecht, the Netherlands. E-mail: vossen@geo.uu.nl 
2 Schlumberger Cambridge Research, High Cross, Cambridge, UK 
3 School of GeoSciences, Edinburgh University, Edinburgh, UK 

Accepted 2004 January 13. Received 2004 January 9; in original form 2003 July 15 

SUMMARY 

Near-receiver material properties are required for the separation of the recorded wavefield into 
upgoing and downgoing P and S waves, and are also important for static time-shift corrections. 
However, it is difficult, especially in land seismics, to obtain reliable estimates for these local 
material properties using conventional techniques. We compare three methods for estimating 
these material properties using a 3-D geophone configuration. The first two methods are based 
on inversion of the wave equation and can be used on almost all of the recorded wavefield. 
However, they require that the wavefield is recorded by a dense 3-D receiver group to allow the 
computation of either spatial wavefield derivatives or interpolants. The third approach is based 
on the inversion of the vertical wavefield propagator. We present a procedure for estimating 
this propagator using only two multicomponent geophones, one buried and one positioned at 
the surface. Propagator estimation and inversion avoids the explicit computation of wavefield 
derivatives, and is therefore less sensitive to measurement errors than both wave-equation 
inversion schemes. However, in the form presented it requires the identification of arrivals 
w> of incoming waves that are isolated in time, and can only be applied to such data. Noise 

tests demonstrate that the propagator inversion provides accurate estimates for P- and S- wave 
velocities of a near-surface low-velocity layer, and is robust with respect to signal-generated 
near-surface reverberations. In case of a near-surface velocity gradient, velocities are obtained 
which are consistent with effective medium velocities. 

Key words: land seismics, multicomponent, near-surface, propagator, wave equation, 
waveform inversion. 



1 INTRODUCTION 

Most observations of seismic waves are made either at, or very near to, the Earth's surface. Before reliable subsurface information can 
be retrieved from these recordings, corrections are required for local near-receiver structure, since variations in this structure often cause 
data perturbations of a similar magnitude to the target signal (e.g. Goupillaud 1961). The variability of the elastic properties close to the 
measurement surface is due to a variety of geological processes and petrophysical properties, among them porosity, permeability, fractures, 
the presence of fluids in pores, compaction, diagenesis and metamorphism (Toksoz et al. 1 976). 

Variations in near-receiver elastic properties cause the following complications. First, receiver static variations in the data are receiver-to- 
receiver traveltime anomalies which occur due to the propagation of most of the seismic energy through the heterogeneous shallow structure. 
Secondly, poor repeatability of the source signature and changes in the source radiation pattern are also attributed to lateral changes in near- 
surface material properties (Aritman 2001). Thirdly, lateral variations in free-surface reflectivity cause differences in the amount of reflected 
and converted energy. This results in amplitude perturbations, especially on horizontal recordings (Kahler & Meissner 1983). Decomposing 
the recorded wavefield into upgoing and downgoing P and S waves allows an analysis of these wavefields without the effects of any free-surface 
interaction (e.g. Dankbaar 1 985; Wapenaar et al. 1 990). However, to perform wavefield decomposition the free-surface reflectivity, and hence 
local subreceiver properties, need to be known. 

Since conventional seismic data are acquired only at the surface, the problem of determining seismic subsurface properties is ill-posed. 
Curtis & Robertsson (2002) therefore proposed to use dense 3-D recording patterns to better constrain land seismic near-surface velocities. 
The pattern consists of a single buried three-component geophone and several surface geophones. The receivers are sufficiently close that 
spatial wavefield derivatives can be computed. These derivatives are required to invert the equation of motion for local material parameters. 
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Robertsson & Muyzert (1999) originally introduced this recording geometry to accomplish P/S separation by explicitly computing the 
divergence and the curl of the wavefield. 

In other research fields, e.g. medical imaging, several concepts to estimate local material properties have been developed. Medical 
practitioners aim to estimate local material properties since variations in the mechanical properties of tissue often reflect early stage pathology 
(Gao et al. 1 996). Unlike in conventional seismic surveys, the displacement field is measured throughout a tissue using measurement techniques 
based on ultrasound or magnetic resonance imaging (MRI) (Muthupillai etai 1995). Quantitative elasticity reconstruction is achieved either 
by comparing modelled stress to measured strain (Gao et al 1996; Van Houten etal. 1999) or by direct inversion of the observed displacement 
field. For example, Romano et al (1998) and Oliphant et al (1999, 2001) have shown that elastic properties can be estimated by localized 
inversion of the equation of motion. 

We propose an alternative method to estimate local near-surface P- and £-wave velocities. This method was originally proposed by 
Trampert et al (1993) to estimate the £-wave velocity structure and the quality factor in a borehole. This was achieved by estimating and 
inverting the vertical SH wavefield propagator derived from the spectral ratio of a downhole data record over a surface data record. We 
formulate this approach for the elastic P-SV case and a strategy for inverting the propagator for near-surface P- and S-wave velocities is 
developed. 

Before discussing propagator estimation and inversion, we briefly review techniques based on direct inversion of the wave equation, 
including those from medical fields. For the purpose of comparison, we illustrate the importance of two of the methods and our propagator 
inversion method using signals perturbed by noise. 



2 WAVE-EQUATION INVERSION 

Oliphant et al (1999, 2001) and Curtis & Robertsson (2002) proposed to constrain local material properties by algebraic inversion of the 
wave equation. A Cartesian coordinate system (*i , jc 2 , * 3 ) is used with the positive vertical direction z = x$ oriented downwards. The particle 
velocity is regarded as a function of space and time t, and is written as v = v(x, t). Overdots are used to indicate time derivatives (e.g. 
v = dv/dt and v = d 2 \/dt 2 ). The wave equation for v in a homogeneous, isotropic medium is then given by (Aki & Richards 2002): 

v = <* 2 V(V ■ v) - p 2 V x (V x v), (1) 

where a = [(k + 2pi)/p] l/2 is the P-wave velocity, ft = {fi/p) Xfl is the S- wave velocity, and X and /x are the Lame parameters. Derivative 
conditions can be derived from the constitutive equation and the free-surface boundary conditions. The constitutive equation for a homogeneous, 
elastic medium reads: 

&ij = kd k v k Sij + M(3fWy + djVi), (2) 
where Sy is the Kronecker delta and 9* denotes the spatial derivative with respect to coordinate x k . Boundary conditions state that the traction 
o ,3 vanishes at the free surface for i = 1, 2, 3, although in practice this is only an approximation, as air waves can exist. By setting a {3 =0 at 
the free surface, we can substitute for vertical wavefield derivative expressions with horizontal derivatives: 

a 3 t>i = -3it;3, (3) 
3 3 „ 2 = -a 2l > 3 , (4) 

v2 . 



(5) 



{a 2 -20 2 \ 

a 3 u 3 = — f — — — i (a ]V| + d 2 v 2 ). 

Curtis & Robertsson (2002) showed that the following system of equations is obtained after substitution of the free-surface derivative conditions 
into the wave equation: 

t), =p 2 A l (t)-a- 2 p 4 B l (t), (6) 
v 2 =p 2 A 2 (t)-a- 2 p 4 B 2 (t) i (7) 

vi=a 2 Ai(t)-p 2 B 3 {t). - (8) 

Expressions for the measurable coefficients A t and £, are given in Appendix A. These coefficients require the evaluation of spatial wavefield 
derivatives. A receiver group as illustrated in Fig. 1 allows the computation of these spatial and temporal derivatives using finite -difference 
operators (assuming that v 2 = 0 and d 2 \ = d 22 \ = 0). The spatial derivatives can only be obtained accurately if the horizontal and vertical 
geophone spacings Ax and Az satisfy the following criteria: Az ~ X™ in /6 and A x ~ AJ lin /6, where A™ n and A™" are the minimum effective 
wavelengths in the x or z direction, respectively (Levander 1988; Muijs et al 2002). Assuming that these spatial derivatives can be calculated, 
eqs (6)-( 8 ) can be inverted for a and p. The attractiveness of this method is that it can deal with the complete wavefield (Curtis & Robertsson 
2002). On the other hand, it requires the computation of second-order derivatives which are likely to be affected by noise and measurement 
errors (Muijs et al 2002). 

Romano et al (1998, 2000) suggested to transform the wave eq. (1) into an integral equation to avoid the computation of spatial wavefield 
derivatives. This approach is referred to as the variational formulation, and we present it for estimating local near-surface material properties 
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Figure 1. Geophone configuration for near-surface velocity estimation. The volume Q bounded by surface r == r b U Tf S with outward-pointing normal vector 
n is used in the variational formulation. Wavefield interpolation is performed along the dashed lines. 



I 



using the geophone configuration shown in Fig. 1 . Consider a vector- valued function w and a volume ft bounded by a surface T with 
outward-pointing normal vector n (Fig. 1). In Appendix B, we demonstrate that eq. (1) is equivalent to 



/ dftw • v = / 



dftw • v 
where 

w Q = a 2 V(V • w) - /3 2 V x (V x w), 



(9) 



(10) 



01) 



w r = (a 2 - £ 2 )[w(V • v) - v(V . w)] - £ 2 [v ■ (Vw) r - w • (Vv) r ]. 

Eq. (9) is known as Betti's theorem (Ben-Menahem & Singh 1981). 

The surface V = T b U r fs , where r fs is the surface of ft which coincides with the free surface (see Fig. 1). Because the vector- valued 
function w can be chosen arbitrarily, boundary conditions can be imposed on w to aid the computation of the integrals in eq. (9). If we assume 
that 



w = 0 

and 
(Vw) T 



on r b 



(12) 



n = 0 on T b , (13) 

the contribution of the surface integral along T b vanishes. The evaluation of the surface integral along r fs requires the computation of 
first-order spatial wavefield derivatives at the free surface. These are also needed for accurate wavefield interpolation throughout ft: first, the 
wavefield and its horizontal derivative are interpolated along the free surface using B-splines (Unser 1999). Secondly, the vector gradient of 
the wavefield is obtained using the free-surface conditions (eqs 3-5). Next, the wavefield in ft is interpolated along lines from the surface to 
the buried geophone (the dashed lines in Fig. 1). The interpolated wavefield is parametrized by second-order polynomials, which are uniquely 
determined by the interpolated wavefield at the free surface, the directional derivative at the free surface and the wavefield at depth Az. For 
sufficiently accurate wavefield interpolation within volume ft, we found that, similar to the previous method, Az ~ X™ n /6 and A* ~ X™ n /6. 
Finally, the surface and volume integrals are evaluated using the trapezium rule. 
Consider the following independent vector- valued functions: 

w^t/^O,*)) 1 *, (14) 



w 2 = (0, 0, f w ) T , 
with 

f w (x h z l ) = {x?-2x? + \)(z l -\) 2 . 

The local coordinates x { and z, are related to x and z according to: 
xi = 2x/Ax n , 



(15) 
(16) 
(17) 



z/ = (2z - Az n )/Az n , (18) 
with Axq and Az n the lengths of the sides of volume ft. The local coordinates are defined on the domain x\ € [—1, 1] and z/ e [— 1 , 1]. The 
functions Wi and w 2 satisfy the boundary conditions (eqs 12 and 13). 

We illustrate both methods with a 2-D example. Synthetic data for a half-space model are computed using the Cagniard-de Hoop method 
(de Hoop & van der Hijden 1983). An explosive line source oriented in the x 2 direction is located at 200 m depth and emits a Ricker wavelet 
with a dominant frequency of 40 Hz. The P-wave velocity is 600 m s" 1 , the S-wave velocity is 200 m s~\ and the density is 1600 kg m -3 . 
Multicomponent geophones are centred around 50 m offset. The horizontal geophone spacing Ax = 1 .0 m and the buried geophone is located 
at 0.50 m depth. For a Ricker wavelet with a 40 Hz dominant frequency, the maximum frequency / max with significant energy is 80 Hz. 



© 2004 RAS, GJI, 157, 796-812 



May 3, 2U04 16:48 Geophysical Journal international gji2249 



^0.395 

CO 




Propagator and wave-equation inversion 799 
(c) (d) 



Figure 2. Synthetic data traces for (a) v x at the surface, (b) v x at depth Az, (c) v z at the surface and (d) v z at depth Az. 
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Figure 3. Sensitivity analysis for direct wave-equation inversion, (a) The waveforms of the left (dashed) and right side (solid line) of eq. (6), and the grey area 
illustrates the region for ±20 per cent perturbations in /J. The dotted line shows the effect of Gaussian noise (25 dB S/N ratio) added to the velocity recordings 
on the right side of eq. (6). Similarly, (b) shows the waveforms of eq. (8) and illustrates the region for ±20 per cent perturbations in ce. The dotted line shows 
the effect of Gaussian noise (25 dB S/N ratio) added to the velocity recordings. 



Then, the minimum wavelength k mhi = c m] Jf m2x = 2.5 m, with c min the minimum phase velocity. The geophone configuration satisfies the 
wavelength condition for the computation of derivatives and interpolants. The length of the sides of Q are Ax a = Ax and Az Q = Az/2 (Fig. 
1). Recorded data are shown in Fig. 2. 

The constraints of direct wave-equation inversion on <x and are illustrated in Fig. 3. For the given model values a and 0, there is a 
good match between the left and right side of eqs (6) and (8). For noise-free data, derivative operators are obtained sufficiently accurately with 
the geophone configuration of Fig. 1. Perturbing a or P results in amplitude changes in the waveforms of the right-side terms in eqs (6)-(8)- 
We compare the sensitivity to perturbations in a and with the effect of measurement errors in the particle velocity recordings. To study 
these effects, 25 dB uncorrelated Gaussian noise (peak-to-peak energy signal-to-noise (S/N) ratio) was added to the velocity recordings. A 
bandpass filter (20 Hz < / < 100 Hz) was applied to these recordings before computing the wavefield derivatives. As a result, the effect 
of noise on temporal wavefield derivatives is reduced. However, spatial wavefield derivatives are still severely distorted by the added noise. 
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Figure 4. Sensitivity analysis for the variational formulation. The waveforms of the left (dashed) and right side (solid line) of eq. (9) are shown for wi (a) and 
w 2 (b). The grey area in (a) illustrates the region for ± 20 per cent perturbations in 0, whereas (b) shows this region for ± 20 per cent perturbations in a. The 
dotted curves show the effect of Gaussian noise (25 dB S/N ratio) added to the velocity recordings. 
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Estimates for a are especially affected, since the wavelength for P waves is larger than for S waves. Consequently, the S/N ratio for spatial 
P-wave derivatives is lower than for S-wave derivatives. 

Fig. 4 shows a similar analysis for the variational approach. There is a good agreement between the left and right side of equation (9) for 
the given model values a and p. In principle, the interpolation and integration steps are sufficiently accurate using the geophone configuration 
shown in Fig. 1 . Perturbing a or p results in amplitude changes in the waveforms of the right side terms in eq. (9). A comparison between 
the sensitivity to a and P and the effects of noise (25 dB S/N ratio) shows that this approach is also very sensitive to the effects of noise, and 
again mostly for a. 

These sensitivity tests demonstrate that it may only be possible to constrain p using these methods. There is not a clear difference between 
the derivative and variational formulation. However, for a geophone configuration with more than one buried geophone, an interpolation scheme 
can be developed which is less sensitive to random errors. 

Muijs et al (2002) studied the effect of deployment-related errors, such as misorientation and mislocation of geophones in a 3-D 
recording geometry, on the computation of divergence and curl. These are particularly sensitive to misorientations of geophones, requiring 
that the orientations of all geophones be accurate within 2°. 

Instead of using the full wave equation, the free-surface condition given in eq. (5) can be used to determine a velocity ratio between a 
and p (Curtis & Robertsson 2002). This only requires the computation of first-order spatial wavefield derivatives, and is therefore expected 
to be more robust than wave-equation inversion. Both the free-surface conditions and wave-equation inversion are applicable to the complete 
wavefield, except the airwave. 

As an alternative, we propose a new method which avoids the explicit computation of spatial wavefield derivatives and imposes no 
requirements on the maximum depth separation of geophones. This is achieved by assuming that the incident wavefield can be described by 
a single plane wave. 



3 PROPAGATOR ESTIMATION 

Our alternative approach to wave-equation inversion is based on propagator inversion. Trampert et al. (1993) proposed estimating the 5-wave 
velocity structure and the quality factor by analysing the vertical SH propagator in the time domain. The propagator contains all information 
on the material properties between the free surface and the depth of a buried geophone, and can be estimated directly from recorded data. 
The S velocity can be constrained by determining the time lag between the peaks in the SH propagator. This time difference represents the 
two-way traveltime of an SH wave to depth Az. The amplitude difference between the peaks can be used to infer the quality factor: downward 
continuation of upgoing waves results in increasing amplitudes, whereas the amplitudes of downgoing waves decrease with depth. 

We generalize this procedure to estimate the elastic P-SV propagator and illustrate a scheme to constrain near-surface P- and S-wave 
velocities based on waveform inversion of this P-SV propagator. 

Assuming a plane wave solution for the wave equation, the plane wave at depth Az can be written in the form 

v(f,x, Az) = P(*,x)*v(*,*,0), (19) 
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where 

(Pn Pn Pn\ 

P=\Pix Pn Pn (2°) 
\Pix Pn Px) 

and * denotes a temporal convolution. For a homogeneous, isotropic and elastic medium between the free surface and depth Az, the propagator 
P is a function of the P-wave velocity a, the S-wave velocity f$ and the horizontal slowness p. For a free-surface incident plane wave with 
slowness p t the theoretical propagator coefficients read (Aki & Richards 2002): 



Ai — B 2 n 2 G p 4- TCI - 26 2 o 2 V21G? 


(21) 


'22 — ^ 1 > 


(22) 


/>„ - rn _ 2£ 2 z? 2 )/21Gf + 0 2 p 2 G? 


(23) 


/> 13 = [/>0 - 2£ V)/(2^)] C7f - £ W<?f , 


(24) 


h\ = P 2 pqpG P 2 - [p(l - 2/* 2 /> 2 )/(2? 5 )] Gf , 


(25) 


Pn = h\ = ^23 = >32 = 0, 


(26) 


where 




Gf(/, />) = [<$(' + «?/>Az) + - Az)], 


(27) 


(/, p) = + tf/>Az) - 5(r - Az)], 


(28) 


GfC P) = + qs Az) + 5(/ - ?5 Az)], 


(29) 


Cf(', />) = P(l + -75 Az) - S(t - q s Az)]. 


(30) 



The vertical slownesses q P and q s are given by: 

q P =(a' 2 -p 2 y^ (31) 

^ = (/r 2 -/> 2 ) 1/2 . (32) 

Eqs (27)-(30) contain the phase shifts for two-way extrapolation of the wavefield towards depth Az. Positive phase shifts describe the 
propagation of downgoing waves to Az, whereas negative phase shifts indicate propagation of upgoing waves to this depth. 

The coefficients of the propagator matrix (eqs 21-26) can be interpreted as follows. Extrapolation of PSV waves is described by Pn, 
Pu, Pi\ and P33 , whereas for SH waves it is given by P n . The amplitude terms before Gf , Gf , Gf and Gf are wavefield decomposition filters: 
the wavefield is separated into P, SV and SH waves before extrapolating the recordings to depth Az, which is given by the phase terms. Finally, 
summation of the extrapolated decomposed wavefield renders the total wavefield at depth Az (Osen et al. 1999; Aki & Richards 2002). 

In the case of a non-attenuating medium, the extrapolation filters for elastic P-SV wave propagation can be obtained directly from the 
data by exploiting symmetry properties of the filters: the theoretical expressions for the extrapolation filters show that P lu P 2 i and F 33 are 
even functions around t = 0 and P\i and Pi\ are odd functions. Hence, the spectra of Pw.Pn and P33 are entirely real, whereas the spectra 
of P| 3 and P31 are purely imaginary. These properties are used to estimate the extrapolation filters without prior information on a, and p. 
Equating real and imaginary parts of eq. (19) in the frequency domain shows that the propagator coefficients are given by: 

Pn = {9*[vj(a>,*. 0)]SR[v,(<w f x, Az)] + 3[v 3 (<», 0)]S[v,(a>, x, Az)]}/D(<w), (33) 

P33 = {M[i>,(a>,x, 0)]9i [113(d), x , Az)] + 3M<y,x, 0)]3[i; 3 (<», x, Az)]}/D(o>), (34) 

P n = i{W[v\(v, *, 0)]S[u,(w,x, Az)] - S[u,(a>,x, 0)]SH[v,(o>, jc, Az)]}/D(o;), (35) 

P31 = *Mu 3 (o>, 0)]3[v 3 (o>, x, Az)] - 3[u 3 (o>, x, 0)Mv 3 (o>, x, Az)]}/0(o>), (36) 
where 

£>(<») = 9t[u 3 (a>,x, 0)Mvi(<w, *, 0)] + 3[u 3 (<w, x, 0)P[v,(w,x, 0)]. (37) 

In these equations, 9t[v(<u, x, z)] denotes the real part of v(co> x, z) and 3[v(<u, x, z)] is the imaginary part of u(o>, x, z). Expressions for P(/, 
x) are found by taking the inverse Fourier transform of eqs (33)-(36). Note that these symmetry properties break down for the P-SV case 
in the presence of attenuation. Then, only the SH case can be treated correctly (Trampert et al. 1993). The propagator estimation procedure 
remains valid in a vertically inhomogeneous medium. 
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The following procedure is used to stabilize the spectral divisions in eqs (33H37). For arbitrary signals F (co) and D(fo), the spectral 
division of F(w) by D(oS) is given by: 
F(<o) 



G(<o) = 



D((o) 



(38) 



Unfortunately, spectral division is numerically unstable both because signals are band-limited and due to the existence of low-amplitude 
notches in the spectrum. In practice, the spectral ratio is estimated using the following technique (Helmberger & Wiggins 1971; Langston 
1979): 



F(o))D*(co) 



(39) 



<t> DD (co) 
where 

<& DD = max{D(co)D*(a)), c max[£>(^)D*(a>)]}, (40) 

and W(a>) is a frequency window to limit the final frequency band in the estimated deconvolution. The complex conjugate of D is denoted by 
D*. The function <b DD can be thought of as simply being the autocorrelation of D(co) with any spectral notches filled to a level depending on 
the parameter c. The frequency windowing function W(co) is a tapered window with cut-off frequencies set to the minimum and maximum 
frequencies for which 

D(co)D*(a>) > c max[D((o)D*(a))l (41) 
This criterion implies that the parameter c also controls the bandwidth of the spectral ratio G'{a)). 
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4.1 Half-space example 

We shall first illustrate the inversion scheme for near-surface velocities with a half-space example. Synthetic data are computed using the 
Cagniard-de Hoop method (de Hoop & van der Hijden 1 983). We restrict ourselves to the P-SV case since SH was fully treated by Trampert 
et al. (1993). The P-wave velocity is 600 m s" 1 , the S-wave velocity is 200 m s -1 and the density is 1600 kg m -3 . An explosive line source 
is located at 200 m depth and emits a 120 Hz Ricker wavelet. Multicomponent geophones are positioned at 50 m offset, one is located at the 
free surface and the second geophone is located at 1 .0 m depth (Fig. 5). Fig. 6 shows traces recorded by these receivers. 

The propagators P u , P\i> Pi\ and ^33 may be estimated from these data. Theoretical solutions for propagator filters are functions of a, 
P and p. Fig. 7(a) shows these filters for the given model parameters. Before comparing these propagator coefficients with the data-estimated 
propagator filters, it is necessary to limit the frequency band of these theoretical expressions to the same frequency band as used for the data- 
estimated propagator filters. The time domain expressions of the frequency windowing functions 0#(a>) (eq. 39), where W i} are the frequency 
windowing functions for propagator components P iJt are shown in Fig. 7(b). Fig. 7(c) illustrates the data-estimated and band-limited theoretical 
propagators. 

We choose to invert the propagator for near-surface velocities in the time domain. In this domain, the propagators can be interpreted as 
follows. The time delay between the peaks in P\ \ and P n gives the two-way vertical traveltime for SV and P waves, and the amplitudes of the 
propagator filters are controlled by the velocity structure between the free surface and depth Az and the signal bandwidth. To constrain a, p 
and p from the propagator waveforms, we define the following objective function: 

E = E U +£33 + £13 + £31, (42) 
with 



50 m offset Az=1.0m 



200 m depth 



5jj a= 600 m s- 1 ; (3=200 m s- 1 ; 

P-source p = 1 600 kg m -3 



Figure 5. Geophone configuration and parameters of the half-space model. 
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Figure 6. Synthetic traces of the particle velocity recorded by the surface and buried geophone shown in Fig. 5. The left graph shows traces of v x and the 
right graph of v z . 



(a) 



Ik 



-o.oi 



0 
t(s) 



0.01 




Figure 7. Propagator fitting procedure for the data displayed in Fig. 6. Given estimates for a, @ and p, theoretical solutions for the wavefield extrapolation 
filters can be computed (a). In practice, data are band-limited, and therefore this theoretical solution cannot be compared directly with the estimated filters, (b) 
The time-domain expressions of frequency filters which limit the solution to the frequency band in which the spectral division was performed. A convolution of 
the theoretical solution in (a) with the filter shown in (b) allows a comparison with the data-estimated propagator. Part (c) shows that there is a good agreement 
between the theoretical (thin lines) and data-estimated propagator filters (thick lines). The stabilization factor c = 1 x 10~ 3 . 

The theoretical solution for the propagator filter component ij is denoted by Py(t t x, a, P, p) and Pij{t, x) is the data-estimated propagator 
filter component ij. The objective function E is a function in a 3-D model space. Cross-sections of the objective function (Fig. 8) show that 
perturbations in p have a relatively small influence on estimates of a and p. 

So far, we have considered noise-free data. The data-estimated propagator filters for recorded data perturbed by 25 dB uncorrelated 
Gaussian noise are shown in Fig. 9. These are compared with the theoretical filters for different values of the stabilization factor c (eq. 40), 
namely c = 10" 4 , c = 10" 3 and c = 10~ 2 . For increasing value of c, the effect of noise is reduced in the propagator filters. Furthermore, the 
frequency band is more limited for a higher c value. The filters are smoothed and contain less energy. 

We estimated the relative variations in a, p and p for different noise levels and c values. To quantify the effect of noise, we define the 
relative root-mean-square error (RMS) in ct by: 



^±(2^2) 



1/2 



(44) 



with a 0 the model P-wave velocity and a, the estimated P-wave velocities, which are found by minimizing E in eq. (42). Similar expressions 
are defined for the relative RMS error in p and p. For each noise level, experiments were conducted 1000 times (N = 1000) with different 
manifestations of Gaussian noise. The minimum of the objective function E was determined using a forward search method. Fig. 10 shows 
that the estimates of a, p and p are most robust for c = 10" 2 : c clearly stabilizes noise very systematically Furthermore, a and p are better 
constrained than p, which implies that relative errors in estimates of p have less influence on estimates of the other parameters than do relative 
errors of similar magnitude in either a or p. 
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Figure 8. Cross-sections intersecting the minimum of the objective function E. Contours are drawn for E = 0.10 and E = 0.20. 
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Figure 9. Effect of the stabilization factor c on noise-perturbed propagator filters. The S/N ratio is 25 dB. The propagator filters are computed for c — 10~ 4 
(left), c = 10" 3 (centre) and c — 10" 2 (right). The thin lines are the band-limited theoretical solutions and the thick lines the data-estimated propagators. The 
amplitudes of each filter are normalized to the solution for c = 10" 4 . For increasing values of c, the noise is effectively reduced. 
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Figure 10. Uncertainties in estimates for or, /J and p for different S/N ratios and c values. 



In addition to random errors, the effect of position and orientation errors of the geophones within the recording pattern is investigated. 
We only consider misorientation in the xz-plane, described by rotation angle 0 (positive for rotation in clockwise direction). The effects of 
misorientation and mislocation are quantified with a relative error e. The relative error in a is given by: 

e(a) = (a - ot 0 )/a Q . (45) 

Similar expressions are used for relative errors in p and p. Fig. 1 1 (a) shows the relative error in estimates for a, ft and p due to misorientation 
of the geophone located at the free surface. The errors in a and P are small compared with relative errors in p. Large errors are found for 9 
close to the angle of incidence. In this example, there are only free-surface incident P waves. This rotation causes the energy on the horizontal 
recording and thus on the denominator D(co) (equation 37) to be minimized. As a consequence, the spectral divisions (eqs 33-36) become 
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Figure 11. Errors in velocity estimates due to mislocation and misoriention: (a) misorientation in the xz plane of a surface geophone; (b) misorientation in 
the xz plane of a buried geophone; (c) vertical mislocation of a buried geophone. The stabilization factor c = 1 0~ 3 . 

unstable. In the case of near-surface structure, we do not expect such a rapidly increasing error, because in that case both P and S waves will 
be recorded. 

The errors caused by misorientation of the buried geophone is shown in Fig. 1 1(b). Here, e{a) and e(p) remain small, because v(a>, Az) 
does not contribute to the denominator D(co). Finally, the effect of vertical mislocation of the buried geophone is illustrated in Fig. 1 1(c) — the 
values for Az are changed in the inversion procedure; the relative change in the depth for the buried geophone Az r = (Az - Az 0 )/ Az 0 , with 
Az 0 the true depth of the buried geophone. There is a linear relationship between errors in Az and estimated velocities. This implies that the 
phase contains most information on the velocities. 

It is difficult to make a fair comparison between wave-equation inversion techniques and propagator inversion. Wave-equation inversion 
uses the amplitudes in the wave equation (eq. 1 ) to constrain a and p. This requires measurement of spatial and temporal wavefield derivatives 
by a dense 3-D geophone configuration, with a spatial separation on a sub wavelength scale. Propagator inversion, on the other hand, uses 
mainly phase information. To estimate time differences reliably, the depth separation between geophones needs to be larger than those used 
for wave-equation inversion. Therefore, we used different depths for the buried geophone and different dominant frequencies for the source 
wavelet. 

The obtained results suggest that propagator inversion is less sensitive to both random and deployment related data perturbations. For a 
S/N ratio of 25 dB, propagator inversion gives no error in a and whereas for the wave-equation inversion techniques Figs. 3 and 4 indicate 
more than 20 per cent deviations in estimates for a, while £ is still well resolved. Furthermore, Fig. 1 1 demonstrated that propagator inversion 
is tolerant to at least 5° misorientations and 5 per cent vertical mislocations of individual geophones. Derivative operators on the other hand 
are particularly sensitive to random misorientations of geophones. Muijs et al (2002) showed that computation of divergence and curl with 
dense 3-D recording geometries requires that the orientations of all geophones are accurate within 2°. This criterion needs to be satisfied as 
well for wave-equation inversion techniques. 

Although propagator inversion is tolerant to measurement errors, it is questionable whether this is the dominant source of errors for this 
method, since a single slowness assumption is implicit in propagator estimation and inversion. Wave-equation inversion techniques do not 
have this assumption and are applicable to the complete wavefield. In the following experiments, we investigate if propagator inversion is 
accurate in a model with near-surface structure. 

4.2 Low-velocity layer example 

The second experiment is performed in a model with a near-surface low-velocity layer. Reverberations in this layer result in multiple arrivals. 
The near-surface layer is 5 m thick with oc = 600 ms _1 ,j6 = 200 m s" 1 and p = 1600 kg m~ 3 . The parameters of the underlying half-space 
are: a = 1500 m s" 1 , p = 400 m s" 1 , and p = 1800 kg m~ 3 . The P-source is located at 100 m depth and emits a 120 Hz Ricker wavelet (see 
Fig. 12). The receiver group is similar to the previous experiment. Synthetic data are computed with a reflectivity method (Kennett 1983). 
The recorded synthetic traces are shown in Fig. 13. These clearly show the multiple arrivals due to interfering waves in the near-surface 
low-velocity layer. 

The data-estimated propagators are again compared with the theoretical propagators (Fig. 14). The latter were computed for the horizontal 
slowness of the first break, which is 3. 1 x 1 0" 4 s m~ 1 . Due to the multiple arrivals (multiple horizontal slownesses), there is not an exact match 
between the theoretical and data-estimated propagators. The energy in the propagator filters decays for higher values of c. For too high a c 
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Figure 12. Model and geophone configuration for experiment with a near-surface low-velocity layer. 




Figure 13. Synthetic traces recorded by the geophone configuration displayed in Fig. 12. Recordings of v x (left) and v z are shown (right). 

value, the fit between the theoretical and data-estimated propagator deteriorates. This can be attributed to the stabilization of internal notches 
in D(co), whereas the theoretical propagator is not compensated for the energy decay in P due to this stabilization. 

Cross-sections of misfit functions using c = 10" 3 are shown in Fig. 15. The minimum is located close to the model velocities, whereas 
the estimated p tends to be smaller than that of the first break, namely around 2.8 x 1 0" 4 sm" 1 . Waves reverberating in the low- velocity layer 
arrive at smaller angles of incidence, thus the estimated p can be regarded as some averaged value over all these arrivals. 

Fig. 16 illustrates the effects of added noise on velocity and slowness estimation. For high S/N ratios, the estimates for a and p do not 
exactly converge for c = 5 x 10~ 3 and c = 10~ 2 , as was suggested by Fig. 1 5 as well. The estimated horizontal slowness differs from that of 
the direct arrival and varies for different values of c. As in the half-space example, c stabilizes the effects of noise. Reliable estimates for a 
and p are obtained for S/N ratios down to approximately 18 dB. 

4.3 Structure between the surface and the buried geophone 

A critical assumption in the propagator inversion scheme is that the medium between the surface and buried geophones is homogeneous. 
However, according to the classical Hertz-Mindlin model (Mindlin 1949), velocities increase with confining pressure and thus with depth, 
especially in the top few metres. This also has been observed in field experiments (e.g. Bachrach et al. 2000). Here, we consider a model with 
a large velocity gradient close to the free surface. The velocity gradient decreases with depth (Fig. 1 7). The gradient model is parametrized by 
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Figure 14. Effect of the stabilization factor c on estimated propagators: c - 10~ 4 (left), c = 10" 3 (centre) and c - \0~ 2 (right). The thin lines are the 
theoretical solutions computed with the horizontal slowness of the first break and the thick solid lines are the propagators estimated from the data. Traces are 
normalized with respect to the c — 1 0~ 4 curves. 
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Figure 15. Cross-sections intersecting the minimum of the objective function E for the low-velocity layer model. Contours are drawn for E = 0.10 and E = 
0.20 and c = 10~ 3 . 
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Figure 16. Uncertainty in estimates for a, £ and p for different c values as a function of the S/N ratio. The estimated horizontal slowness is compared to the 
slowness of the direct arrival. 



0.02 m thick layers. Synthetic data are computed using a reflectivity method (Kennett 1983). The explosive line source is located at a depth 
of 1 00 m and 50 m offset, and emits a 120 Hz Ricker wavelet. The buried geophone is located at a depth of 1 .0 m. Fig. 1 8 shows the synthetic 
traces for both the surface and buried multicomponent geophones. Data-estimated propagators are shown in Fig. 19. These are compared 
with theoretical solutions with effective medium velocities, the lower and upper bounds of which are given by the Reuss and Voigt average 
velocities (Wang & Nur 1992). For the top 1.0 m, the Reuss velocities are a" = 259 m s _1 and p~ = 101 m s" 1 , and the Voigt velocities are 
a+ = 284 m s" 1 and f} + = 105 m s~ ! respectively. The average value of the Voigt and Reuss velocities is used to compute the theoretical 
propagator. Although the propagator waveforms do not exactly match, the time lags of the peaks are similar. Fig. 20 shows cross-sections of 
the misfit function. Both a and ft are well constrained, whereas p cannot be resolved. The estimated velocities are a = 272 ± 8 m s" 1 and 
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Figure 17. Gradient velocity model and source location. 
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Figure 18. Synthetic data recorded in the gradient model at 50 m offset. The buried geophone is located at a depth of 1.0 m and the P source at a depth of 
100 m. 

p = 102.8 ± 0.6 m s _l . These values are obtained by performing repeated experiments (N = 1000) with different manifestations of random 
noise (S/N = 25 dB), determining the minimum of the objective function E with a forward search method, and using the average and standard 
deviation of all these estimates. The obtained velocities fall within the Voigt and Reuss bounds for effective medium velocities. 



5 DISCUSSION 

Three methods to estimate near-surface properties were evaluated in this paper. The first two methods are based on inversion of the wave 
equation, whereas the third method inverts wavefield propagator filters for near-surface velocities. A fundamental difference between the 
methods is that wave-equation inversion schemes use amplitude information to constrain a and By changing a and the amplitudes in 
the wave equation change, not the phase. On the other hand, propagator inversion uses predominantly phase information to constrain a and 
/J, i.e. the time lags of the peaks in the propagator coefficients are controlled by the wave velocities. 
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Figure 19. Data-estimated and theoretical propagators in the gradient model. The theoretical propagator is computed using the average value of the Voigt and 
Reuss effective medium velocities: a = 272 m s" 1 and 0 = 103 m s" 1 , p = 8.75 x 10 -4 (s m" 1 ) and c - 10" 3 . 
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Figure 20. Cross -sections through the minimum of the objective function E for the gradient model. A contour is drawn for E = 0.20 and c = 10 3 . 



These sensitivities for either predominantly amplitude information or phase information result in a different optimum experimental 
set-up. Wave-equation inversion techniques require the measurement of spatial wavefield derivatives, with a spatial geophone spacing on a 
subwavelength scale. Propagator inversion on the other hand uses mainly phase information. To estimate time differences accurately, both the 
depth separation of the surface and the buried geophone and the bandwidth of the recorded signal have to be larger. For this reason, we used 
different configurations, different frequencies and a different depth of the buried geophone, to compare wave-equation inversion techniques 
with propagator inversion. 

Before fitting the waveforms in the wave-equation inversion schemes, either spatial wavefield derivatives or interpolants need to be 
evaluated. For an accurate implementation of these methods, the spacing between geophones must be sufficiently close, approximately 1 /6 
of the effective wavelength. Since the wavelength depends on the material properties, some prior knowledge has to be available to design the 
receiver group. 

Propagator inversion on the other hand has no wavelength constraints on the maximum allowed separation of geophones. Therefore, 
the method may be applicable to borehole recordings. This has been demonstrated for the SH case by Trampert et al. (1993). However, in 
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the P-SV case we only treated the elastic situation, whereas attenuation has to be included in the formulation for borehole applications. The 
near-surface velocity estimates presented in this paper are not affected by attenuation, because the geophone separation Az <3C A.. 

However, for a geophone buried at too shallow a depth, propagator inversion cannot resolve the wave velocities. Then, small errors in 
the phase will cause large uncertainties in the estimated velocities, especially for ct. These uncertainties are reduced by increasing the depth 
of the buried geophone (Van den Berg et al. 2003). The minimum depth for the buried geophone depends on the frequency band, the angle of 
incidence, on a and on the time sampling rate. 

The wave-equation inversion schemes use the complete wavefield as input signal. Any event which satisfies the wave equation provides 
constraints on the near-surface properties. The main sources of error which need to be dealt with are deployment-related errors. Muijs 
et al. (2002) demonstrated that errors such as misorientation and mislocation of individual geophones significantly affect estimates of spatial 
wavefield derivatives. Correction schemes need to be developed to compensate for these effects. Propagator inversion, on the other hand, 
avoids the explicit computation of spatial wavefield derivatives. As a consequence, this method is less sensitive to these types of errors. 

In contrast to wave-equation inversion techniques, a single slowness assumption is required to construct and invert the propagator 
filters. As a consequence, those data containing arrivals that are isolated in time need to be selected before this method can be applied. 
We demonstrated, however, that in a medium with a near-surface low-velocity layer, estimates for a and /J are not significantly affected by 
signal -generated reverberations, and that in a model with a strong velocity gradient, effective medium parameters are obtained. 

The following strategies for propagator inversion can be developed by combining data of different shots. First, if the horizontal slownesses 
in selected data show little variation for different shot points it is possible to reduce effects of noise by stacking the estimated propagator 
filters. Second, propagator inversion can be formulated without requiring the isolation of an event which can be approximated by a single 
plane wave. When using a dense source array, the frequency-wavenumber spectra can be computed for common receiver gathers of both the 
surface geophone and the geophone at depth. Consequently, propagator filters can be determined for each wavenumber or horizontal slowness 
individually by a spectral division of these two spectra, Thus, propagator inversion is potentially applicable to the complete wavefield. However, 
difficulties might arise due to lateral variations in near-surface material parameters and poor repeatability of the source. 

Synthetic tests have been performed in two dimensions, whereas 3-D wave propagation has to be considered in the field. Applications 
in the field require a different acquisition geometry for wave-equation inversion. Wave-equation inversion techniques require a 2-D patch 
of 4 x 4 receivers at the surface and a buried geophone in the centre at this geophone group. Then, all spatial and temporal derivatives or 
interpolants can be computed. Propagator inversion requires the separation of P-SV waves from SH waves, hence rotation of the recordings 
in the in-line/cross-line direction is required. Furthermore, propagator inversion is aided by an estimate for the horizontal slowness in the 
in-line direction, and therefore a cross pattern of receivers at the surface is recommended for field applications. The general characteristics of 
the considered methods do not change in three dimensions, and therefore the presented 2-D synthetic results will apply to the 3-D situation. 

6 CONCLUSIONS 

We evaluated three methods for estimating local near-receiver material properties. A fundamental difference between the methods is that 
the wave-equation inversion schemes constrain the /'-wave and S-wave velocities by matching amplitudes in the wave equation, whereas 
^^^^ propagator inversion constrains the wave velocities using predominantly phase information, i.e. with the traveltime from the free surface to 
fly li t the depth of the buried geophone. As a consequence, wave-equation and propagator inversion have different optimum receiver configurations. 
I II Si fll ^ or wave-equati 00 inversion, the configuration needs to be designed to allow the measurement of spatial wavefield derivatives or interpolants, 
OBiJ whereas the configuration in propagator inversion must allow a reliable measurement of a phase difference. 

Both the derivative and the integral formulation for wave-equation inversion are almost equally sensitive to the effects of measurement 
errors. If a configuration with more than one buried geophone would have been allowed, an interpolation scheme may be developed which 
is less sensitive to these type of errors. However, these configurations are not convenient for practical applications, and therefore were not 
considered. 

Propagator inversion has two advantages and one disadvantage over wave-equation inversion schemes. First, instead of a 3-D receiver 
configuration, only two geophones are required, one positioned at the surface and one buried. Second, this scheme avoids the explicit evaluation 
of spatial wavefield derivatives or interpolants and therefore it is less sensitive to deployment-related errors. On the other hand, this technique 
implicitly assumes that a data window containing arrivals with a similar horizontal slowness can be isolated in the recorded data, whereas wave- 
equation inversion schemes are applicable to nearly complete recordings. Synthetic experiments were performed to assess the consequences 
of the single slowness assumption in propagator inversion. These experiments demonstrated that the method is robust with respect to signal- 
generated reverberations, and in the case of a near-surface velocity gradient, results are consistent with effective medium velocities, the upper 
and lower bounds of which are given by the Voigt and Reuss averaged velocities. Moreover, when using a shot array, propagator inversion 
could be formulated without the plane wave assumption, and therefore it is potentially applicable to the complete wavefield. 
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APPENDIX A: COEFFICIENTS FOR DIRECT WAVE-EQUATION INVERSION 



The measurable coefficients in eqs (6)-{8) for direct wave-equation inversion are found by inserting the free-surface derivative conditions, 
eqs (3}-{5), into the wave eq. (1): 



2 

A { (t) = —(9^3 + hv\) + V£ u, + 23, (V„ ■ v„), 


(Al) 


^(0 = -^"(32 "3 + 33 K2) + V> 2 + 23 2 (V„ • V„), 

Az 


(A2) 


M0 = ^(V w • v w + Bjuj) - V 2 „v 3 , 


(A3) 


5,(0 = 23, (V„ ■ v„), 


(A4) 


S 2 (/) = 23 2 (V„ ■ v„), 


(A5) 


and 




fi3(0=^(V w v w )-2(V^ 3 )- 


(A6) 



The finite-difference (FD) first-order derivatives in depth are denoted 3 3 v and are given by: 



3 3 v(Az/2)= V(Az) - V(0) + O(A^), 
Az 
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with Az the depth beneath the free surface of the buried geophone. The vertical derivative cannot be measured exactly at the free surface but 
at depth Az/2. The following notation is used for horizontal wavefield derivatives: 

v* = [a, a 2 ] T , (A8) 
and 

y H = [u, v 2 ] T (A9) 
Temporal and spatial wavefield derivatives are computed using finite-difference operators. 

APPENDIX B: INTEGRAL EQUATION FORMULATION 

We demonstrate that the wave eq. (1) can be transformed into an integral equation. Let w = (w \ , w 2 , wi) T be a sufficient smooth vector- valued 
function and £2 a volume bounded by a surface T. If the wave equation is multiplied by w and integrated over the volume the result is: 

f yv\dQ = f w [a 2 V(V • v) - /3 2 V x (V x v)]dQ. (Bl) 
Jn Jn 

The first term of the right-hand side of eq. (Bl) can be expanded by repeated application of the divergence theorem: 

/ dQvi • [V(V . v)] = / <*r[w(V • v)] • n - f dQ(V • w)(V ■ v) 
Jn Jr Jn 

= j dT[w(W • v) - v(V ■ w)] ♦ n + J dQv ■ [V(V • w)], (B2) 
where n is the outward-pointing normal vector on V. The second term is expanded using the identity, 

V x (V x v) = V(V . v) - V 2 v (B3) 
and the divergence theorem: 

j dQw . [V x (V x v)] = J dr[vt(V . v) - v(V • w) - w • (Vv) T + v ■ (Vw) T ] • n 

+ / ^^v • [V x (V x w)]dQ. (B4) 

In the volume integrals, all spatial wavefield derivatives have disappeared in eq. (B4). Because the vector-valued functions w can be chosen 
arbitrarily, boundary conditions can be imposed such that the surface integrals vanish. Consequently, all spatial wavefield derivatives can be 
made to disappear. 
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